Metric tensor as the dynamical variable for variable cell-shape molecular dynamics 
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We propose a new variable cell-shape molecular dynamics 
algorithm where the dynamical variables associated with the 
cell are the six independent dot products between the vectors 
defining the cell instead of the nine cartesian components of 
those vectors. Our choice of the metric tensor as the dynami- 
cal variable automatically eliminates the cell orientation from 
the dynamics. Furthermore, choosing for the cell kinetic en- 
ergy a simple scalar that is quadratic in the time derivatives of 
the metric tensor, makes the dynamics invariant with respect 
to the choice of the simulation cell edges. Choosing the densi- 
tary character of that scalar allows us to have a dynamics that 
obeys the virial theorem. We derive the equations of motion 
for the two conditions of constant external pressure and con- 
stant thermodynamic tension. We also show that using the 
metric as variable is convenient for structural optimization 
under those two conditions. We use simulations for Ar with 
Lennard- Jones parameters and for Si with forces and stresses 
calculated from first-principles of density functional theory to 
illustrate the applications of the method. 

PACS: 61.50.Ah, 71.15.Pd, 62.20-x, 62.50+p 



I. INTRODUCTION 

With the development of new simulation methods and 
the increase in available computational power, molecular 
dynamics has become an important tool in the simulation 
of matter in the condensed state |[],^| . In its earliest ap- 
plications, molecular dynamics methods were employed 
to simulate systems of interacting particles with a con- 
stant density and energy, using a simulation cell with a 
fixed volume and shape and with a constant number of 
particles inside. For extended systems, periodic bound- 
ary conditions were introduced to reduce finite cell-size 
effects. 

The calculations with constant energy, volume and 
number of particles are expected to simulate the ther- 
modynamic properties of the micro-canonical ensemble. 
However, in laboratory conditions, one often controls the 
intensive variables temperature T and pressure p, instead 
of the extensive variables E and V. Therefore molecular 
dynamics methods were developed to simulate systems 
at constant temperature or pressure fl|-[j]]. In the case 
of constant pressure simulations, the size and shape of 



the simulation cell must be allowed to change. In order 
to do so, an "extended system" is constructed which in- 
cludes degrees of freedom for the cell. A microscopic sim- 
ulation of the structural, mechanical, and dynamical re- 
sponse of material systems to external stress of interest in 
tribology, material fatigue and wear, crack propagation, 
stress induced phase and structural transformations, lu- 
brication and hidrodynamical phenomena, is more con- 
veniently done with varying cell shapes. 

The dynamics of the cell is fictitious. Therefore there 
are many reasonable choices for the equations of motion 
of the variables associated with the cell. Traditionally 
the dynamical variables were the cartesian components of 
the vectors defining the periodicity of the simulation cell. 
The early choices for the equations of motion had some 
invariance problems, and more complicated equations of 
motion have been proposed to avoid those problems. 

Here we suggest the use of the dot products between 
the vectors defining the simulation cell as the variables 
for the cell dynamics. We show that using the new vari- 
ables avoids in a natural way the problems previously 
encountered. 



II. VARIABLE CELL-SHAPE MOLECULAR 
DYNAMICS 

To simulate a system at constant pressure, one must 
allow for variations of the volume and shape of the sim- 
ulation cell. Andersen || proposed to use the volume V 
of a cubic simulation cell as a dynamical variable in an 
extended hamiltonian, thus allowing for volume fluctu- 
ations driven by the dynamical imbalance between the 
imposed external pressure, p cx t, and the actual instanta- 
neous internal pressure, p- m t, as given by the virial theo- 
rem. As the simulation cell is periodically repeated, the 
dynamics associated with the cell is fictitious. In the ex- 
tended lagrangian for the dynamics, Andersen included 
a fictitious kinetic energy term associated with the rate 
of change of volume, 
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where W A is a fictitious "mass" associated with the cell. 
He also added the term U ce \i = p ey ±V , which is the poten- 
tial from which the constant external pressure acting on 
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the cell is derived. During the simulations, the volume V 
fluctuates about an average value such that, in the limit 
of long simulation times, the time average of the calcu- 
lated internal pressure is equal to the chosen external 
pressure, p int = p ey ±. Here we use an overline to indicate 
the limit of a time average for long calculation times. In 
those simulations it is the enthalpy H = E + pV that 
is approximately conserved, not the internal energy, and 
Andersen showed that, assuming ergodicity, his simula- 
tion method samples the isoshape-isobaric-isoenthalpic 
ensemble to an accuracy of O (N~ 2 ) when calculating 
ensemble averages of intensive parameters {O (A -1 ) for 
extensive parameters), where N is the number of parti- 
cles in the simulation cell. 

Andersen's method is best suited to study equilibrium 
properties of fluids, for which the shape of the cell is 
irrelevant. To study shear flow (viscosity) in fluids or 
to study solids it is not enough to change volume with 
constant shape. For example, a given cell shape may be 
compatible with the periodicity of one crystal structure 
and be incompatible with another solid phase, and so the 
fixed cell shape may artificially prevent the appearance of 
thermodynamically more stable phases. In order to study 
structural phase transitions, Parrinello and Rahman |i]|7]] 
extended Andersen's method to allow for changes in both 
the volume and the shape of the cell. They used as dy- 
namical variables the cartesian components 

of the three vectors Sj defining the periodicity of the 
simulation cell. Here efj are the three orthonormal vectors 
that define a cartesian coordinate system. To generate 
the dynamics, a fictitious kinetic energy of the cell 

W PR 3 3 

i=i j=i 

is included in the lagrangian, where W PR is again a fic- 
titious mass. In the limit of large N, the equipartition 
principle tells us that the kinetic energy of the 9 variables 
of the cell is small compared with the kinetic energy of 
the 3iV — 3 variables associated with the particles' posi- 
tions, and the method simulates the isobaric-isoenthalpic 
ensemble. 

As the kinetic energy of the cell is fictitious, it can be 
chosen in many reasonable ways that simulate the same 
ensemble in the limit of large number of particles, N, and 
large simulation times. However, different choices of the 
fictitious cell kinetic energy yield different dynamics, and 
one can ask which is better or more convenient. Several 
authors have pointed out some shortcomings of the origi- 
nal method of Parrinello and Rahman: it is not invariant 
under modular transformations (defined below), the con- 
sistency between the condition of mechanical equilibrium 
and the virial theorem is only verified in the large N limit, 
and it has spurious cell rotations ff-|ll]]. 



For a given periodic system, there are infinite equiva- 
lent choices of the basic simulation cell. If a, are three 
vectors commensurate with the periodic system, then the 
transformation a'^ = ^2 k Mkjtik, with M an integer ma- 
trix with |detM| = 1, gives another set of vectors de- 
scribing the periodicity. It is desirable that the dynamics 
should not depend on the particular choice that is made, 
i.e., the equations of motion should be formally invari- 
ant with respect to the interchange between equivalent 
cells (modular transformations) H,||. This characteristic 
improves the physical content of the simulation, by elim- 
inating symmetry breaking effects associated with the 
fictitious part of the dynamics ||. Of course, in the ther- 
modynamic limit (N — > oo) these effects vanish, but they 
may be important in computer simulations, which may 
use only a small number of particles. That is often the 
case in first-principles molecular dynamics (l2| . 

For long simulation times and constant applied pres- 
sure, the dynamics for the cell should yield (Peart)} = 

Pext$j, with (Peart)} the internal stress in cartesian co- 
ordinates given by the stress theorem, which is a gener- 
alization of the virial theorem [p^[ . A weaker condition 
that is easily checked is that this should be verified in 
particular when the cell is restricted to undergo isoshape 
fluctuations p0| . Andersen's method obeys this condi- 
tion, while the same is true for the Parrincllo-Rahman 
dynamics only in the large N limit Q . 

The orientation in space of the simulation cell is irrel- 
evant for the structural and thermodynamical descrip- 
tion of the system (principle of material-frame indiffer- 
ence [fil)). However, it is included in the dynamics if one 
uses the components of the cell edges as dynamical vari- 
ables, and spurious cell rotations have been obtained in 
actual simulations with the Parrinello- Rahman method, 
namely in the simulation of molecules, whose internal de- 
grees of freedom sometimes cause the internal stress to 
be asymmetrical jl4| . These rotations not only are phys- 
ically irrelevant, but may complicate the analysis of the 
simulations' results. Methods to eliminate them have 
been proposed, such as constraining the matrix of the 
lattice vectors to be symmetrical ]14[ or upper triangu- 
lar Jig] (geometrical constraints), or by symmetrization 
of the infinitesimal strain at each time step (dynamical 
constraint) [ jTl| . 

III. USING THE METRIC AS A DYNAMICAL 
VARIABLE 

If Si, a,2 and c?3 are three linearly independent vectors 
that define the periodic simulation cell and form a right- 
handed triad, then all the properties of the simulated 
system depend only on the symmetrical metric tensor, 

and not on the orientation of the three vectors in space. 
In our new simulation method, we use the 6 independent 
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components of the metric tensor as the dynamical vari- 
ables for the cell. The three diagonal elements of the 
metric give information about the lengths of the lattice 
vectors, and the three independent off-diagonal elements 
contain the additional information about the angles be- 
tween those vectors. The covariant components of the 
tensor g are related to the matrix h = (01,02,03), the 
transformation matrix between cartesian and lattice co- 
ordinates, by the relation 

9 = h T h, (2) 

where h T is the transpose of h. The one-forms b l associ- 
ated with the lattice vectors Si, which are (except for a 
factor of 2ir) the reciprocal lattice vectors, are related to 
the contravariant components of the metric tensor, 

9 v ee // • V = 

and the volume of the unit cell is given by 

V = det h = ^/detgij. 

The position r(i) of the i th atom in the simulation cell 
can be defined by its lattice coordinates s k (i), 

f(i) = s k (i)a k , 

where we use the Einstein summation convention for ten- 
sorial quantities. The distance between any two points 
can be calculated from the metric and the lattice coordi- 
nates, and therefore they completely define the geometry 
of the simulation cell. 

In the Parrinello- Rahman formalism, the calculation of 
the total distance traveled by a particle can be mislead- 
ing. Because the unphysical motion due to the rigid rota- 
tion of the cell should be discarded, one cannot in general 
simply use a two-point formula to calculate, for example, 
the mean square displacement of a particle. The correct 
formula for the distance is naturally expressed in terms 
of the metric: 

As = / J s l g i:j s3dt. (3) 
J t 

For a fixed cell (i.e., for a fixed gij), the Newton's equa- 
tions of motion can be derived from the lagrangian 

C 1 (s i (k),s i (k),g ij ) = 
k^2k=i m ( k )^( k )9i^ j (k)-U(s i (k),g ij ), [ 1 

where the summation is over all N atoms in the cell, 
m{k) is the mass of atom k, and the potential energy per 
cell U includes interactions between atoms in different 
cells. U is a function of the 3N lattice coordinates of 
the atomic positions and the 6 independent components 
of the metric tensor. In the examples of a latter sec- 
tion, U is either the Born-Oppenheimer energy from a 
first-principles pseudopotential local-density calculation 



or the potential energy of the Lennard-Jones model. The 
momentum canonically conjugate to s I (k) is 

= g.f^ = m(% ij s J (fc), 

and the corresponding hamiltonian is 

Wx (a<(fc), = £ Vi{k 2^ {k) + U ('*(*).«*) ■ 

(5) 

To construct the extended lagrangian for the cell dy- 
namics, we must choose the fictitious kinetic energy of the 
cell, -Keen, and, for simulations with applied pressure, add 
the term p ey ±V = p e xt J det g[j . A simple non- negative 
scalar that is quadratic in the time derivatives of all the 
components of g is 

K° cc ii(9v,9^) = — [f t ) [f t ) = —9n {9 lk 9kl9 b ) 

where W° is a fictitious cell "mass" which has the dimen- 
sions of mass times length squared. The positivity of this 
term is shown in Appendix A. Instead of K® cll , we choose 
the slightly modified expression, which is again a scalar, 
quadratic in g, but with a different densitary character, 

^c S di i9ij,9ij) = ~Y (^9ij)9ji {g %k 9kig l3 ) , 

where W s is a fictitious cell "mass" with the dimen- 
sions of mass times length -4 . Alternatively, we may 
view M.i lkl = W 9 (det gy) g lk g 1 ^ as an effective "mass" 
tensor. Although K^ ell gives slightly more complicated 
equations of motion for the cell, it has the advantage of 
reducing to Andersen's K^ n (see Eq. 0) in the case of 
isoshape fluctuations of the cell, if we make the identifica- 
tion W 6 = jW A , and so the dynamics that it generates 
obeys the virial theorem in that limit |Icj . Since we are 
using the metric, the orientation of the cell never appears 
in the equations. It can also be verified that A^ cll is in- 
variant with respect to modular transformations of the 
a l . 

The fictitious lagrangian for the extended system in 
the presence of an applied external pressure is 

£2 {s i {k),s i {k),g i j,gi^) = \ J2 k m(k)s l (k)g lj s 1 (k)- 
U (s l (k),gij) + ^ {<te,g i j)g ji g ik g k ig l: > - p cxtv /det g~ h 

(6) 

and the equations of motion for the atomic coordinates 
are 

m(k)s l (k) = gVFjW ~ m(k)g» 9]l s l (k), (7) 

where Fj(k) = —dU/ds^k) are the covariant compo- 
nents of the force (which can be viewed as the compo- 
nents in reciprocal lattice coordinates multiplied by 2w) . 



3 



This equation for the scaled atomic coordinates is iden- 
tical to the one obtained from Parrinello-Rahman's la- 
grangian, since it does not depend on the choice of K ce \\. 
It should be stressed that this doesn't imply that the 
dynamics of the atoms is the same, because in order to 
convert from the scaled dynamics to the actual atomic 
dynamics we have to use the metric, which is determined 
by the cell's dynamics. Hence the importance of a fic- 
titious cell dynamics which doesn't introduce unphysical 
symmetry-breaking effects || . 

The coupling of the atomic motion to the cell's motion 
is made through the second term on the R.H.S. of Eq. 
[?], which is independent of the orientation and state of 
rotation of the cell; from this, the physical irrelevance of 
the orientation of the cell is evident. 

The equation of motion for the cell variables is derived 
with the help of the relation -^^Aetgij = g lk detgij, 



giving 



dot Qi j \ ^/ dct Qi 

{g lk g kl gij - g kl g k ig tj ) + ^ {g k ig lm g mn g nk ) g tj , 

(8) 

where the contravariant components of the internal stress 
are (see Appendix B) 



V ij = V m(k)sHk)sHk) ~ 2^—. 
k dgi i 



(9) 



The instantaneous internal pressure averaged over the 
cell is gijTr'Pj, and it can also be obtained from Tii 
or C\. 



Pint = 



fdHi 



dCi 



Defining the momentum canonically conjugate to the 
metric tensor 

BC 

n« = — 1 = w* (det^) <fW j = n* 

ogij 

the conserved extended hamiltonian can be written as 



H 2 (s i (k),g ij ,TT i (k),U^)=j: k - 



U(s l (k),g i3 ) 



2W% dct g; 



i(fc)7r'(fc) 
2m(k) 

PextV. 



+ 



(10) 



IV. ANISOTROPIC EXTERNAL STRESS 

A constant applied anisotropic stress is in general non- 
conservative, and thus there is no conserved extended 
hamiltonian in a constant anisotropic stress simulation 
P,^6[. Of course some experimental situations are es- 
sentially non-conservative, and therefore best simulated 
by an appropriate non-conservative dynamics p 16 1. In 



this section we will present a conservative dynamics, but 
one should keep in mind that the simulation should be 
taylored to the problem. 

Molecular dynamics simulations with an applied 
anisotropic stress were first proposed by Parrinello and 
Rahman [Q] . Ray and Rahman |l7| later showed that the 
original formulation was valid only in the limit of small 
deformations, and they proposed an extension valid for 
finite deformations, in which it is the thermodynamic ten- 
sion (defined below), not the stress, that is kept constant, 
and the quantity that is approximately conserved during 
the simulation is the generalized enthalpy of Thurston 
]lc| . This approach is based on the fact that, if the ex- 
ternal stress is allowed to change when the cell deforms, 
so as to keep the thermodynamic tension constant, the 
virtual work of the stress upon deformations of the cell 
is conservative, and so that stress is derivable from a 
potential, which can be used to construct an extended 
hamiltonian. The thermodynamic tension is given by 



V 

To 



hoh 



-1 cart 



'0 1 



(ii) 



where ho and Vo are the reference lattice and its volume, 
and (Text* i s the external stress in cartesian coordinates. 
For h = ho, t and cr^t* coincide. The virtual work SW 
done by an external stress on the faces of the cell during 
an infinitesimal deformation of the cell in the state h is 



5W = V Tr (rde) 



(12) 



where e is the strain tensor for the lattice h measured 
from the reference lattice ho- We see that Vqt is the 
thermodynamic variable conjugate to the strain. Thus, 
for fixed t, the differential is exact, and so we can inte- 
grate 8W over a finite deformation, to obtain the elastic 
energy 



In the following sections it will be convenient to define 
a symmetrical covariant internal stress tensor as 



U cen (h)= / 5W = V Tr(Te). 

Jho 



int " %i' 

which contains the contributions from the potential en- 
ergy UtoV. 



The generalized enthalpy of Thurston is given by |T^] 

H = E + V Q Tr (re) , 

where E is the energy of the system. For our metric- 
based formulation, it is desirable to use the metric, in- 
stead of the strain, as the thermodynamic variable. In 
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order to find what is the conjugate variable, we have 
to express SW in terms of infinitesimal variations of the 
metric tensor. This can be done for a symmetrical (i.e., 
torque-free) external stress, which does no work in pure 
rotations of the cell. The result is given in Eq. (3.5) of 
Rcf |l7jl , and, expressed in tensorial notation is a simple 
expression, 



SW = -Tr 
2 



cxt<5 9 jk 



The thermodynamic variable conjugate to the metric is 
therefore the external stress in contravariant lattice co- 



ordinates. Keeping a, 



*J7 
ext. 



oilt constant when the cell 



deforms thus leads to a conservative external stress, de- 
rived from the potential 



Ucen (g) 



1 



5o, 



(13) 



where go = hyho is some reference metric, 
fixed, one can drop the constant term 
the definition of U ce \\, obtaining 

^ceU {g) = ^<4lt9ij, 



Since <T^ t 



1 % 3 

2 CT cxt50i 



IS 

from 



(14) 



which is independent of a reference configuration and 
quite compact when compared with the definition of r 
in Eq. |ll]. 

The condition that er^t is constant is equivalent to re- 
quiring t to be constant because its cartesian coordinates 
are 



Tab 



1 , 



as can be seen using Eq. ^ from Appendix B jl9[ . Nev- 
ertheless, to a given a 1 ^ doesn't correspond a unique 
thermodynamic tension, because ho is arbitrary. All the 
physical information is contained in c^t an d gij-, ex- 
cept for the (arbitrary) choice of axes. The thermody- 
namic tension fixes the choice of axes and also a refer- 
ence state, through ho- Notice that from the transforma- 
tion law for the contravariant components of the stress, 



dct h 
dct h 



ha^th 1 



obtained from Eq. 



PQ keeping tr^S* constant, we can conclude that U cc i\ as 
given by Eq. U4 is invariant under modular transforma- 
tions ||. 

To grasp the physical meaning of er^ti let's consider 
the force acting on the face of the cell opposing edge i. 
Using Eqs. |lj and ||, we obtain (in cartesian coordinates) 



r = - 



dU, 



cell 



(15) 



showing how the force on the face i is related to the stress. 

The new extended lagrangian can be obtained from 
C 2 given by Eq. |, by replacing p cxt V by the new U cc i h 



Eq. |TJ. The equation of motion for the atoms, Eq. [?], 
remains unchanged, and the equation of motion for the 
cell is obtained from Eq. || by replacing p e xt<7y by -^of**: 



W£g l3 



2 dct gij 



w* 



(gikg kl gij 



— [gkig g mn g ) g%j- 



The conserved hamiltonian is 



iiso 

U{s i (k),g ij ) + 11 



2ffs dct gi 



EiTi(k)iz z (k) 
k 2m(fe) 
I 1 

' 2 a cxtgij- 



l gkigij) 



(16) 



(17) 



In specific applications, it may be desirable to impose 
a constant external pressure, p cy ±, plus a constant ther- 
modynamic tension. Note that the stress tensor associ- 
ated with a constant pressure is a l e J xt = PcxtVg 1 - 1 , and 
so constant pressure is not a particular case of constant 
thermodynamic tension. The generalization is straight- 
forward, and in this case, when considering only isoshape 
fluctuations of the cell, the equation of motion for the cell 
becomes 



V 



where Wa is Andersen's cell mass. This equation shows 
that the off-diagonal elements of (V — er C xt)j are re- 
stricted to be zero, and the diagonal elements are re- 
stricted to take equal values, at all times: by imposing 
a fixed cell shape, we have arrived at an isotropic total 
stress, as should be expected on physical grounds. In 
equilibrium V = 0, and so the average of each diago- 
nal component of yVj equals y (o" cxt )* + p C xt$j, which 

implies = p~~^ + jyTr (a ext )], where the R.H.S. is 
the total external pressure. This shows that our method 
obeys the virial theorem in the case of isoshape fluctu- 
ations of the cell (the proof in Ref. |10| mentioned in 



Section III was for applied pressure only) 



V. STRUCTURAL OPTIMIZATION 

A problem encountered in the simulation of materials is 
the determination of the equilibrium structure of a crys- 
tal at a given pressure (or anisotropic stress) predicted 
by a given model U (s l (k),gij) of its total energy. This 
can, in principle, be achieved by the minimization (under 
the appropriate constraint) of U, which is quite difficult 
because it is a multivalleyed function of many variables. 
A practical strategy is to use a simulated annealing to 
bring the configuration to a deep valley, followed by a 
search of a minimum in that valley. The annealing step 
can be carried out by the variable cell shape molecular 
dynamics described previously coupled to a thermostat, 



5 



brownian dynamics forces, or a periodic rescaling of the 
velocities. The local minimization can be done efficiently 
if one has the gradient of the function to be minimized. 

If we want to obtain the crystal structure at zero tem- 
perature and for an applied pressure of p ex t> we must 
minimize its enthalpy, 

H (s l (k),gij) = U (s l (k),gij) +p cxt y / det g li . 

The gradient of the enthalpy with respect to atomic po- 
sitions is 



dH 



dU 



ds^k) ds l {k) 



which is minus the covariant components of the force on 
that atom. Notice that in molecular dynamics it is the 
contravariant components, F l (k) = g^Fj(k) that appear 
in the equation of motion. The gradient of the enthalpy 
with respect to the metric is 



dH 



dU 







.. = t. l-p ex t~ — \J det . 

<>!/.. <>g,j og lj 



1 



1 



- o a int + 77Pcxt.9 y v^Sy 



The minimum is obtained when the forces are zero and 
when the mixed stress tensor divided by the volume is 
the pressure times the identity tensor, 



f 



\/det g~j 



as desired. 

If we want to obtain the crystal structure for a fixed 
thermodynamic tension, then we must minimize the gen- 
eralized enthalpy, 



f 



H {s l (k) 7 g l3 ) = U (s\k) l9t3 ) + -a e i t g ir 



(18) 



The gradient of the generalized enthalpy with respect to 
the atomic lattice coordinates is still minus the covariant 
force on the atoms, and the gradient with respect to the 
metric is 



dH 



dU 
dgtj 



1 



1 



1 



which is zero when the internal stress is equal to the 
desired applied stress, a\^ t = o"ext- 



VI. APPLICATIONS 

In this section, we apply the method to the study of 
structural phase transitions and structural optimization. 
The isobaric-isoenthalpic ensemble, besides being some- 
what unusual ^0|, is not the most adequate to study 
transitions induced by pressure or stress, because it does 
not allow for the exchange of heat with the surroundings. 
There are several methods described in the literature to 



perform simulations at constant temperature by connect- 
ing the system to a "heat bath" (|. !]. In our examples we 
use Langevin molecular dynamics \ 4], but first we checked 
that in our simulations the generalized enthalpy was con- 
served in the absence of the heat bath. In the course of 
a simulation of a structural transformation the release 
of heat of transformation must be dissipated to the heat 
bath, this takes some time and therefore the temperature 
of the system may rise to values quite above those of the 
heat bath. 

As our first example we simulated a silicon crystal un- 
der a constant pressure using first-principles molecular 
dynamics. Pioneering examples of first-principles molec- 
ular dynamics with variable cell shape include the opti- 
mization of the structure parameters of MgSiC>3 under 
pressure pl| and the structural transition of silicon un- 
der pressure p2|| . In the first case the dynamical variable 
for the cell was the strain tensor, in the second the lattice 
vectors. 

In our simulation of Si, the energy, forces and stresses 
were calculated within the local-density approximation, 
'using a pseudopotenti al p 3fl and a plane- wave basis set 
with a cutoff of 16 Ry pifTThe simulation cell contained 
8 atoms, initially disposed in a diamond structure, with 
lattice constant a = 9.435 a.u. The applied pressure was 
25 GPa. The equations of motion were integrated with a 
Beeman algorithm pj|. The time step was h = 200 a.u., 
and the cell "mass" W 9 = 10 a.u. Langevin dynamics 
with a viscosity damping constant of 7 = m ^ Si ^ , where 

m (Si) is the atomic mass of silicon, was used to simulate 
a heat bath with a temperature of 300K. 

It is well-known that silicon undergoes several 
phase transformations with increasing pressure, and its 
pressure-volume phase diagram has been extensively 
studied (^6|. Starting from a diamond lattice, the struc- 
ture changes at ~11 GPa into /3-Sn, and between 13 
and 16 GPa transforms into simple hexagonal. Other 
densely packed phases appear at around 38 GPa. In the 
first ~ 0.7 ps (200 steps) of the simulation, we observed 
(Fig. 1) that the volume of the simulation cell was fluc- 
tuating around a value that corresponds to the volume 
of the metastable diamond structure of Si at 25 GPa 
(V ~ 885 a.u. for the 8 atoms of the conventional cubic 
unit cell). There was then a rapid drop in the volume, 
accompanied by a rapid rise in the ionic temperature 
to around 3500K (well above the melting point). The 
simulation was interrupted after 1000 steps, well before 
equilibrium with the thermal bath was reached. After 
the transition, the volume of the cell oscillated around 
650 a.u., slightly below the volume of the stable simple 
hexagonal structure at that pressure, but above the den- 
sity of the close-packed structures. Remembering that 
at atmospheric pressure Si contracts upon melting, and 
considering the high temperatures of the simulation, our 
results indicate that at high pressures, the liquid phase 
may still be denser than the solid phase. 

For the purpose of illustrating a molecular dynamics 



6 



method the origin of the forces is irrelevant, therefore we 
used a Lennard-Jones model, with the constants adjusted 
to simulate argon, for the other examples in this article, 
as the computational demands are much lower. In our 
second example, we started with a cubic simulation cell 
with 32 argon atoms in an fee lattice, and increased one 
diagonal contravariant component of the external stress, 
cr^t, linearly with time from zero to 15 x 10~ 5 a.u. for 
the first 4000 simulation steps, and held thereafter cr^ t at 
that value. All other applied stress components were kept 
at zero. This corresponds to a situation of uniaxial com- 
pression. During the first 10000 steps of the simulation 
the system is kept in contact with a heat bath at 10 K, at 
which point we minimize the generalized enthalpy (Eq. 
ILsf ) using a method by Davidon |2j| . The minimization 
is obtained in 96 steps which is approximately the num- 
ber of variables, indicating that the heat bath kept the 
system near a quadratic region of the potential. 

During the compression the system yields for an ap- 
plied stress of ~ 0.1 GPa (after ~ 2500 steps), and due to 
the rearrangement of the atoms, the applied stress drops 
to a minimum of ~ 0.07 GPa and then rises again grad- 
ually, as the thermodynamic tension is increased. After 
the structural rearrangement, the argon is still in a dis- 
torted fee lattice, but the stress is now applied in a [110] 
direction instead of the initial [100] direction, and the 
area on which the force is applied is ~ \[2 times larger. 
The yield was accompanied by a rapid rise in the ionic 
temperature up to ~ 33 K. The heat was gradually dis- 
sipated, and at around step 4000 the temperature was 
back to 10 K. 

Figure || shows the evolution of two of the contravari- 
ant lattice components of the internal stress, cr^t an d 
CTg^t i compared with the corresponding imposed external 
stress components. At first the internal stress oscillates 
around the external values, in particular it accompanies 
the rise in applied stress. When the system yields we ob- 
serve a dramatic increase in the amplitude of the stress 
oscillations, which are then damped with time. Finally, 
in the minimization step the internal and external stress 
are identical within the precision demanded in the mini- 
mization (10 -5 ). 

The contravariant components of the stress tensor are 
not what we are used to call stress (their dimensionality 
is energy per area), so we show in Figure |^ the evolution 
of the cartesian components er zz of the applied and inter- 
nal stress, where we chose the z axis to be in the direction 
of the applied stress. One can see that the cartesian com- 
ponents of the applied stress are not constant when the 
contravariant components of the stress are constant, and 
that the oscillations of the internal stress are magnified, 
but they track each other, and they are identical at the 
end of the enthalpy minimization, as desired. 

The yield is also apparent in the plot of the poten- 
tial component of the generalized enthalpy (Eq. |l8|) as a 
function of time (Fig. |4|). First we observe an increase 
of the enthalpy during load due to the work done on the 



system by the uniaxial stress. When the system yields 
there is a strong decrease of the potential component of 
the enthalpy even while we continue loading the system, 
showing that energy is transfered to the kinetic compo- 
nents, and later dissipated to the heat bath. Only near 
the end of the loading cycle do we see the enthalpy rising 
again. During the annealing steps there is a rapid initial 
decrease of the enthalpy, meaning that the minimization 
procedure rapidly reaches the valley of the multi- variable 
function, but then takes some time to reach the mini- 
mum. 

The best way to observe the yield is from the plot of 
the evolution of the lattice constants (Fig. ^|). From the 
initial slopes one could extract the Young and Poisson 
moduli for the system. After the yield we see that the 
three lattice constants are different from each other, and 
that they are rapidly determined by the minimization 
procedure. At the end of the simulation and after the in- 
spection of the angles we obtain a monoclinic simulation 
cell, which is in reality a supercell of the orthorhombic 
system one should expect when loading an fee crystal in 
the [110] direction. 

A movie of the simulation shows that the {100} planes 
parallel to the uniaxial stress become distorted close- 
packed planes by compression along the direction of the 
stress and expansion along the perpendicular direction. 
A similar simulation was performed by Ray and Rahman 
in Ref. |p8| . They found an fee to close-packed transition, 
with the final structure presenting stacking- faults. 

Our final example is of a structural optimization under 
pressure. We start from conditions quite away from equi- 
librium, perform 2000 steps in contact with an heat bath, 
and then switch to a gradient minimization. Our target 
pressure is 0.3 GPa and we simulate 16 argon atoms with 
a Lennard-Jones potential. The final structure is close- 
packed and corresponds to a stacking of close-packed 
planes that is neither fee nor hep. During part of the 
simulation the temperature is well above melting, so the 
memory of the initial configuration is lost. The evolution 
of the potential contribution to the enthalpy is shown in 
Fig. [| In the inset of that figure that magnifies the 
minimization part of the simulation, one can see that we 
obtain the enthalpy of Lennard-Jones argon at that pres- 
sure. The true minimal structure is not reached because 
the energy cost of the stacking faults is too small, so the 
procedure only finds a deep local minimum. 

In principle the calculation of the energies, forces and 
stresses can be carried out within the metric formal- 
ism, and therefore one never needs to construct the lat- 
tice vectors, that is the matrix h. Our code for the 
Lennard-Jones interaction was written to test the present 
formalism and is fully implemented in the metric lan- 
guage. It never uses the matrix h. Our pseudopoten- 
tial plane-wave code is based on Sverre Froyen's Berke- 
ley code, which stored atomic positions and fc-vectors 
in lattice and reciprocal lattice coordinates respectively, 
that is in contravariant and covariant coordinates. In 
fact the stress was calculated by applying the chain rule 
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dU/dhij = {dU / dgki)(dgki / dhij) , so it was easy to con- 
vert the program to the present formalism. The calcu- 
lations involving the separable non-local pseudopotential 
projectors are easier to perform in cartesian coordinates, 
so for that specific case we construct from the metric g 
a triangular h and proceed in cartesian coordinates. The 
arbitrary choice of the orientation of h has, of course, no 
effect in the results of the calculation. Our plane-wave 
code also has a old, but convenient, symmetry recognition 
package that only works for the conventional orientation 
of the unit cell. If one wants to perform simulations with 
fixed symmetry, than one has to put "by hand" the de- 
sired orientation of h before using the package. Replacing 
those two parts of the code to avoid using the matrix h 
is a straightforward, but tedious job, it is much easier to 
use the tested old subroutines and construct a matrix h 
whenever it is needed. 



VII. CONCLUSIONS 
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IX. APPENDIX A 

To prove that -K^eii * s non-negative, we have to show 
that gjig lk gkig b = Tr (gg^gg^ 1 ) is non-negative. Using 
Eq. |^ and the usual properties of the trace, we find that 

Tr {gg- 1 gg- 1 )=2Tr[(X + X T )X], 

where X = hhr 1 . Writing the rightmost X as 
\ (X + X T ) + \ (X - X T ) and using Tr (X T X) > 0, 
we arrive at the desired result. It will be useful to derive 
this result in the hamiltonian formalism: 

Defining Hba.ij = ^It ^ J ; where some h compatible 
with g was chosen, and defining 



We have shown that the metric is a very convenient 
dynamical variable to use in molecular dynamics simu- 
lations with variable cell-shape. As the cell part of the 
dynamics is fictitious, there is no unique choice of the 
kinetic energy, to be included in a lagrangian or hamil- 
tonian formulation. The use of the tensorial notation 
in a metric formalism, with the requirement that the 
energy functions must be scalars, restrict our choice of 
those functions. The simplest expression for the cell ki- 
netic energy has several properties that were not present 
in early expressions, namely absence of rigid rotations 
and invariance with respect to modular transformations. 
With a convenient choice of the densitary character of 
the kinetic energy, the virial theorem is also satisfied for 
isoshapc fluctuations. For anisotropic stress, the simplic- 
ity of Eq. [l4] contrasts with the definition of thermody- 
namic tension, (Eq. [ll]) and its dependence on a reference 
cell. 

From our kinetic and potential functions for the cell 
metric, we derived the equations of motion for variable 
cell shape molecular dynamics under the conditions of 
constant applied pressure and anisotropic applied ther- 
modynamic stress. We also showed that the optimiza- 
tion of structures under both conditions can be naturally 
expressed in the metric language. Simulations of sili- 
con with first-principles forces and argon with empirical 
Lennard-Jones forces were used to illustrate the applica- 
tions of our equations of motion and minimization pro- 
cedures to the study of systems under applied pressure 
or stress. 
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where ab = H a b,kh we can write the kinetic energy of 



the cell as 
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where P a b = H a b,kl^ = Pba is a new generalized mo- 
mentum for the cell. The canonically conjugate coordi- 
nate is Q ab = \Hab,ij) 9ij> as can be seen using the 
Poisson brackets relations between canonically conjugate 
variables (in this last expression, H T is viewed as a 9x9 
matrix with indices ab and ij). The relation between the 
variables (Q,P) and the variables (<7,n) is similar to the 
relation between ("7*, 7?) and 

From Eq. [II], it is clear that K 9 ell is positive and con- 
tributes with 6 distinct quadratic terms to the energy, 
and so the equipartition theorem applies to the degrees 
of freedom of the cell when they are in contact with a 
heat bath. 



X. APPENDIX B 

The stress is not a true tensor, but a tensorial density 
p9[ , thus transforming differently from tensors under a 
change of coordinates whose jacobian is not unity. The 
transformation of a tensorial density T> lJ from cartesian 
to lattice coordinates is given by |29] 
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V kl = (deth) h^V, 



(20) 



The average symmetrized internal stress in cartesian 
coordinates is obtained from the stress theorem [|13l : 



8 



' cart 



dE 



X; o m(a)v l (o))J 3 (a) 



=o 

dU 



(21) 



where E is the internal energy, v l (a) is the velocity of the 
atom a, and e is the (symmetrical) Lagrangian strain 
corresponding to a rotation-free infinitesimal homoge- 
neous deformation given by h = 1 1 + e ) h J30| , from 



the state g to the state g . In order to convert to lat- 
tice coordinates, we use Eq. |2^ and apply the chain rule 

= ir^ together with the relation g - g = 2h T e h 



0. The result is 



V 13 = -2 



dE 
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; dK 

' dg i3 



, dU 



(22) 



where K and U are the kinetic and the potential energy of 
the atoms in the cell, respectively, and the factors of two 
arise from formally treating gij and gji as independent 
variables. The kinetic energy is given by the first term 
in the R.H.S. of Eqs. || or |5|. Using these expressions, we 
find, with the help of the relation 
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and so it must be made clear whether it is the {iTi (k)} 
or the {s l (k)} that are kept constant when taking the 
derivatives in Eq. It is easily seen that, in order 

to obtain the correct kinetic internal stress (see Eq. |2l| ) 
when transforming Eq. 22 back to cartesian coordinates, 



the {TTi(k)} must be kept constant, and thus we conclude 
that 



-pi] 
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dgi. 



= 2 



fdd 



y / 7r m (fc) 
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which is Eq. ^. 

We note that, because 



\J dct 9ij 



\ d 9ij. 
is a scalar capacity j29| , 



and Vij is a tensorial density, the product , Vjj 

that appears in Eq. |g is a true tensor, and thus has 
the same densitary character as the pressure term in the 
same equation, p ox tgij- 
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FIG. 1. The volume (in atomic units) of an eight-atom 
Si cell is shown as a function of the step of a first-principles 
molecular dynamics with an applied pressure of 25 GPa. The 
dots indicate the simulated data and the three horizontal lines 
indicate the volumes of the diamond, simple hexagonal and 
face centered cubic structures of eight Si atoms at a 25 GPa 
pressure. At that pressure the stable phase is the eight-fold 
coordinated simple hexagonal. The volume starts by oscil- 
lating around the volume of the initial diamond phase, but 
after 200 steps shows a rapid decrease to values near the 
equilibrium value at 25 GPa. The heat released during that 
transformation melts the system and at the end of the short 
simulation it has not yet reached equilibrium with the sur- 
rounding heat bath. The simulation seems to indicate that at 
a pressure of 25 GPa, like at zero pressure, Si contracts upon 
melting. 

FIG. 2. Two of the contravariant (lattice) components 
of the applied and internal stress tensors (in atomic units) 
are shown for a simulation of a cell of 32 argon atoms with 
a Lennard- Jones pair potential submitted to uniaxial load- 
ing. One of the diagonal components of the applied stress is 
increased from zero to 15 x 10~ 5 a.u. during the first 4000 
simulation steps and held constant thereafter, while all the 
other components are held at zero. During the first 10000 
steps the system is kept in contact with a heath bath at 10 K. 
Thereafter we minimize the generalized enthalpy as described 
in the text. As the minimization is very fast, the horizontal 
scale is multiplied by a factor of 20 in that region. During the 
molecular dynamics the internal stress (wiggly lines) oscillate 
around the applied stress (straight lines) as it should. The 
minimization makes the internal stress equal to the external 
stress within the tolerance of the minimization procedure. At 
~ 2500 molecular dynamics steps the system yields in the way 
described in the text. 

FIG. 3. The cartesian components of the applied and 
internal stress tensor along the direction of compression are 
shown for the same simulation as in Fig. 2. The applied stress, 
which is only indirectly controlled through its contravariant 
lattice components, also oscillates. In particular, during the 
phase transformation the applied stress drops considerably in 
response to the decrease in the average internal stress due to 
the atomic rearrangement. 

FIG. 4. The potential component of the generalized en- 
thalpy (in atomic units) is shown as a function of time for the 
same simulation as in Fig. 2. First we observe an increase of 
the enthalpy during load due to the work done on the sys- 
tem by the uniaxial stress. When the system yields there is 
a strong decrease of the potential component of the enthalpy, 
even while we continue loading the system, and only near the 
end of the loading cycle (indicated by the arrow) do we see the 
enthalpy rising again. The horizontal scale is again multiplied 
by a factor of 20 in the minimization part of the simulation, 
showing the efficiency of the procedure of enthalpy minimiza- 
tion in an expanded scale. 



FIG. 5. The evolution of the three lattice constants (in 
atomic units) with time is shown for the same simulation of 
the previous figures (Figs. 2, 3 and 4). The strong structural 
rearrangement during yielding is clearly seen. 

FIG. 6. The enthalpy (potential part only and in atomic 
units) of a cell with 16 atoms of Lennard- Jones argon is shown 
for a simulation with an applied pressure of 0.3 GPa. The 
simulation starts at conditions quite away from equilibrium, 
evolves for 2000 steps in contact with a heat bath, and then 
the enthalpy is minimized. The inset shows the minimization 
part of the simulation in an expanded scale. The horizontal 
line in the inset is the enthalpy of Lennard- Jones argon at 
0.3 GPa, and that is the value reached by the minimization 
procedure. Dots that seem out of place in the minimization 
correspond to overshooting steps in the multi-dimensional 
minimization procedure. 



10 




MD Step 



Fig. 1 Souza 




-0.0001 1 1 1 1 1 1 1 1 1 1 1 1 

2000 4000 6000 8000 10000 10100 

MD Step 

Fig. 2 Souza 




MD Step 



Fig. 3 Souza 



-0.08 1 1 1 1 1 1 ■ 1 1 r 




2000 4000 6000 8000 10000 10100 

MD Step 

Fig. 4 Souza 




2000 4000 6000 8000 10000 10100 



MD Step 



Fig. 5 Souza 




Fig. 6 Souza 



